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Abstract. Old maps are often the basic source of modern environmental 
studies focused on landscape reconstruction. Georeferencing of these maps 
is a key activity in the implementation of these studies. It is necessary to 
know the theory of the basic methods of georeferencing. This paper 
summarizes the current knowledge of methods, and further proposes a new 
method for the georeferencing of multiple sheet map series. It is based on 
the overall adjustment of all transformation coefficients with conditions of 
the map edges continuity. 
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1. Introduction 

Recently, old maps are increasingly used as a base layer for various 
environmental applications. With the expansion of map servers, web map 
services and APIs of mapping portals they are increasingly used as an es- 
sential base layer for monitoring changes in the landscape or as a source of 
information for the reconstruction of the original landscape. For a long 
time, the biggest problem with the inclusion of old maps in these 
appi i cati ons i s bei ng sol ved; their georef erenci ng. 

At the department of mapping and cartography we are working on the 
problematics of georeferencing of old maps for a long time. Some results of 
our research can be found in (Krejci 2009) or (Cajthaml 2010). 

1.1. Global Transformation Methods 

There are many methods of georeferencing of old maps. Core group of these 
transformations are global transformation methods. They are based on a 
unique transformation key applied to the entire map area. The coefficients 
of these transformations are usually determined by adjustment by elements 



using the least squares method. The advantage of this method is to deter- 
mine the most probable values of the coefficients including their 
characteristics of accuracy. According to the transformation equations we 
can find the following most commonly used methods: linear 2D 
transformations, polynomial transformations. 

2D linear transformations are based on linear relations between the local 
coordinates on the old map and global coordinates in a defined coordinate 
system (identical points). In principle, these relations consist of set of 
simple geometric operations, resulting in a corresponding transformation. 
A map can be moved, rotated or its scale can be changed. Of course, these 
operations can be carried out independently in both coordinate axes. That 
results in the following transformation methods: similarity, 5-parameter 
affine, affine, and projective. The relevant transformation equations can be 
found in most works dealing with georeferencing. (Beineke2001, Livieratos 
2006) 

Polynomial transformations, unlike the linear transformations, use a poly- 
nomial relation between the coordinates. Generally, in practice, the second 
and third order of polynomials is used (quadratic and cubic). Higher poly- 
nomials do not significantly improve the transformation method. 

Besides the transformations that are defined by equations for the entire 
region, there is the possibility to transform the coordinates using methods, 
where each point corresponds to a separate mathematical relation. These 
methods can be classified as local, because the transformation key changes 
throughout the enti re area. 

1.2. Local Transformation Methods 

The principle of local transformation methods are relationships based on 
the position of the transformed point. Distance to identical points is very 
often used as a parameter. This is actually the analogy to the interpolation 
in the space between theidentical points. Identical points then usuallyshow 
zero residuals after the transformation. Transformations are non-residual, 
unlike the previous ones, where just residuals of identical points are used 
for least squares adj ustment. 

Global methods preserve the spatial relationships in the image by the trans- 
formation equations (depending on the method used) and local methods 
distort the image so that the identical points were non-residual. There is 
essential to cover the enti re area with identical points with appropriate den- 
sity for local methods. Especially areas outside the envelope of identical 
points can be much distorted. On the other hand for the old (incorrect) map 



local transformations can better fit the image to the current state (e.g. for 
layers overlap). 

Local transformation methods are widely used in combination with global 
methods. In the first step global transformation coefficients are computed; 
in the second step residuals after the transformation are interpolated. Two 
of the best known and most widely used methods include: inverse distance 
weighted (IDW) (Shepard 1968), thin pi ate spline (TPS) (Bookstein 1989). 

Local transformation methods can also include transformations that are 
calculated by the i mage segments. The most common way is to divide space 
into irregular triangles whose vertices consist of identical points. There is a 
certain analogy with modeling surfaces using TIN. If the space is divided 
into triangles, there is need to select the appropriate transformation model. 
The simplest and most common method is to use the linear model. In each 
triangle there are determined six unknown parameters of affine transfor- 
mation using coordinates of three apex points. The advantage of this proce- 
dure is primarily in its simplicity. Properties of transformed areas corre- 
spond with affine transformation. At the intersection of two triangles the 
map drawing is continuous (property of affine transformation), but lines 
are not smooth. This is due to the fact that the applied transformation sur- 
face is not smooth (itisgeneral polyhedron). Besides using the linear model 
it is possible to transform triangular areas with smooth surfaces. In this 
case, we require the continuity of the transformation surfaces of individual 
triangles. There is necessary to ensure that the surfaces have always com- 
mon tangent plane at the apex. The theory of these surfaces is beyond the 
scope of commonly used methods. More can be found in (Farin & Hansford 
1977). 

Another option is the division of the map by the quadrangles. Such nets are 
not much used in practice; however, they can have advantages. This method 
can be used i n cases where we have i denti cal poi nts on the map frame. Then 
it is advantageous to divide the entire image by these frame marks. Each 
part of the image is then transformed separately and uneven map defor- 
mation is removed. Use of such method is to remove the map sheets' 
shrinkage. Projective transformation is used as a transformation method. It 
fits precisely using four identical points (corners). Use of projective trans- 
formation is the disadvantage of this method, where map drawing at the 
intersection of two quadrangles after transformation can be discontinuous. 
Most of map deformations are not so big that the discontinuity is not no- 
ticeable, but sometimes it can occur. Then it is necessary to look for other 
solutions using local non-residual transformations. 



2. Map series georeferencing 

When georeferencing maps, we can easily encounter situations where we 
are georeferencing two or more adjoining map sheets. Then it is desirable 
that the edges of adjacent map sheets exactly match. There are several pos- 
sibilities how to ensure the continuity. The basic method is to merge the 
whole image before georeferencing. It can be used with a small number of 
map sheets and relatively good quality of image data. Using graphical soft- 
ware the image data can be combined so that they fit exactly on the edges. 
The large amount of resulting data (practically unrealistic is to merge a 
larger number of map sheets) is limiting this method. Other problems can 
be need of good knowledge of graphic software for the quality result or the 
inability to simply merge data (map sheets) with different coordinate sys- 
tems. 

Another option has two steps. I n the first step the coordinates of the cor- 
ners of map sheets are determined and in the second step the projective 
transformation using these corner coordinates is performed. Projective 
transf ormati on i s advantageous i n terms of transformi ng j ust four i denti cal 
points (corners) when it becomes non-residual. First, it is necessary to de- 
termine the coordinates of the corners of map sheets. This can be done ei- 
ther directly (from the knowledge of map sheet frame) or indirectly from 
identical points. Separate transformation of each map sheet can get the dif- 
ferent coordi nates of the same corner poi nts. One corner poi nt mi ght corre- 
spond to a total of four adjacent map data sheets, and thus may also be de- 
termined four times. Then, the resulting coordinates are determined as an 
average. The averaged coordinates are then possible to be used within pro- 
j ecti ve transf or mati on . 




Figure L Map sheets after separate affine transformations (left) and after 
following averaging of corners and projective transformations (right). 



As mentioned, the use of projective transformation is the disadvantage of 
this method, where map drawing at the intersection of two quadrangles 
after transformation can be discontinuous. Two step transformation, which 
can be time consuming, can be another problem. It is also not possible to 
transform the map with curved edges (only a linear shape of the frame is 
possible). Relatively easy procedure without a deeper knowledge of graphic 
software or a more complex adjustment could be the main advantage. 
Another method that will be described in more detail and presents newly 
proposed solution is overall adjustment of transformation coefficients of all 
map sheets with the conditions that define the edges continuity. The meth- 
od works with identical points of all map sheets and with conditions defin- 
ing edges connectivity. Here a calculation of all the coefficients is performed 
simultaneously. The image data are then transformed separately according 
to determined coefficients. This eliminates the dual transformation of im- 
age data and worki ng with huge data. Unfortunately, there is need for a new 
adjustment when changing any identical point, because everything is com- 
puted together. 




Figure 2. Map sheets after separate affine transformations (left) and after 
new proposed method with overall adjustment (right). 



A more complex method based on surfaces geometry is using parametric 
surface patches. It is possible to define an area of patch using boundary 
curves (map sheet edges). It is necessary to obtain identical points on these 
edges. Then it is possible to reconstruct a new patch shape. This method is 
applicable to maps with frame markers, through which you can get good 
shape of boundary curves. (Coons 1977) 



2.1. Affine transformation with edges continuity 

As a f i rst exampl e of overal I adj ustment of transf ormati on coeff i ci ents there 
are i ntroduced relations for the affine transformation, which isprobablythe 
most widely used method. For simplicity, the selected base case consists of 
the transformation of two adjacent map sheets. Left map sheet has number 
land the other (right) map sheet has number 2. 

For sheet No. 1 we have transformation: 

x' - aix + hy + X tl , 
y' = c x x + d x y + Y tl . 

For sheet No. 2 we have transformation: 

x' = a 2 x + b 2 y + X t2 , 
y' = c 2 x + d 2 y + 1<2 ■ 

Conditions of identical shape of both common edges are based on the iden- 
tity of two common corner poi nts (ends of the edge). Thus: 
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where the coefficients \1 are map sheets numbers, indices L, R represent 
the left, respectively right corner and indices T, B are top, respectively 
bottom corners of the map sheet. Substituting the transformation equations 
we get: 



o i -' nr + b^yRT + X n = a 2 2 xlt + b 2 y LT + X t2 , 
CiXrt + d^ynr + Y n = c 2 2 x LT + d 2 2 y LT + Yii , 

a/avw + b\ l y RH + X n = a 2 2 x ul + b 2 y LH + X t2 , 
c^xrb + dxy RB + Y tl = c 2 2 x LB + d 2 y LB + Y t2 . 



Condition equations after small change are: 

a^XRT + b L l y RT + X n - a 2 2 x LT - b- 2 2 y LT - X t2 = . 

c^Xrt + di l y RT + Y a - c 2 2 x LT - d 2 2 y LT - Y t2 = . 
ayXfui + bi l tjRB - A',.i - a 2 x U} - b 2 2 y LU - X l2 = . 

ciXrb + di-ijRB + Y tl - c 2 x LB - d. 2 2 y lb - Y t2 = 



Now classic method of adjustment by elements with constraints can be 
used: 
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M atrix B consists of condition equations derivatives: 
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I i s the vector composed of def i ned coordi nates of the i denti cal poi nts i n the 
target system. After a simple calculation we can determine all unknown 
transformation parameters for the two map sheets. Condition equations 
ensure that thecorners fit perfectly after the transformation. 

For the successful adjustment it is necessary to define only independent 
conditions. When adjusting multiple map sheets there must be previously 
think about which conditions should be used. In the particular case, for four 
map sheets touching at one point it is possible to use only 3 conditions of 
equality coordinates from a total of 6 possible combinations of pairs of co- 
ordinates (map sheets). The remaining conditions can no longer be used, 
because they are linear combination of three used conditions. 

Transformation coefficients applied to the images of the map sheets ensure 
the identity of both shared corners of map sheets and map content is abso- 
lutely continuous. If we wanted eliminate distortion of the map frame (e.g. 
bended linear frame), it would be necessary to reconstruct the image shape 
before georeferenci ng. This can be achieved by the aforementioned method 
of surface patches. 

If we want to use proposed method with a different (easier) 2D linear trans- 
formation instead of affine, one should bear in mind the geometric nature 
of this method. This means that before the computation all map sheets 
should be able to fit together after used transformations. 



2.2. Polynomial transformation with edges continuity 

A second example of overall adjustment of transformation coefficients can 
be the second order polynomial transformation. It was chosen mainly be- 
cause of relatively frequent use of this method for georeferenci ng of old 
maps. Partly similar method proposes Molnar (2010). Unfortunately, his 
method works only with perfect rectangular maps, cropped by the bounding 
box parallel with coordinate axes. In fact, his method is special (simpler) 
case of our new general method. 



For simplicity, a base case will be shown again (transformation of two adja- 
cent map sheets). Left map sheet has number land the other (right) map 
sheet has number 2. 



For sheet No. 1 we have transformation: 

x' = (ijx 2 + Ihy 2 + cixy + diX + e x y + f\ , 
y' = 9\x 2 + hjy' 2 + iixy + j x x + fay + h . 



For sheet No. 2 we have transformation: 

x' = a 2 x 2 + b 2 y 2 + c 2 xy + d 2 x + e 2 y + f 2 , 

y' = ff2^ 2 + ^y 2 + i2%y + + k 2 y + ^ ■ 



I n the case of using a polynomial transformation the identity of the corner 
points is not enough for ensuring continuity through the whole common 
edge. The general solution of the common edge identity is the identity of 
functional form of the edge curve. I n practice, the identity of three points of 
the curve (for second polynomial order) after transformation is sufficient. 
Conditions of identical shape of both common edges are: 
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where the coefficients \2 are map sheets numbers, indices L, R represent 
the left, respectively right edge and indices T,B,M are top, bottom, 
respectively middle points of the edge. Substituting the transformation 
equations we get 6 condition equations. 



M ethod of adj ustment has the same scheme as the affi ne method. I n addi- 
tion, matrices A and B have similar shape with polynomial transformation 
derivatives for A and condition equations derivatives for B. Of course the 
size of matrices is different due to different number of transformation pa- 
rameters. 

After calculation we can determine all unknown transformation parameters 
for the two map sheets. Condition equations ensure that all the edges fit 
perfectly after the transformation. Transformation coefficients applied to 
the i mages of the map sheets ensure that the map content i s absol utely con- 
tinuous. 



3. Experiments and discussion 

As the most suitable candidate for the use of the newly proposed method 
seem to be results of the first military mapping survey of Austria-Hungary 
(1763-1785). There were selected two testing areas where identical points 
were collected. After their statistical testing they entered into the calcula- 
tion of affine transformation with edges continuity. For both areas the 
results are very promising. All map sheets are continuous and fit precisely. 
Detail of four adjoining sheets can be seen on the Figure 3. 

The second method with using second order polynomial transformations 
will be tested in the near future. Another plan to the future is evaluating of 
proposed georeferenci ng methods precision. When using least squares 
method it is possible to compute standard errors for all determined param- 
eters as well as for anyidentical point. Using this evaluation it is possible to 
describe quality of the georeferenced map sheets. First tests show that 
standard positional error of the whole produced map is two times bigger 
than standard positional error of the individual map sheet. Comparing both 
methods (affine and second order polynomials) show that using polynomial 
method reduces the positional error by half. 

There seems to be a very interesting project of georeferenci ng of all map 
sheets of the first military mapping survey in the area of Czech Republic. It 
also seems possible to include proposed algorithms in the software 
Georeferencer (Pridal 2011) for map series georeferenci ng. 



Figure 3. Detail of four neighboring map sheets after using the new pro- 
posed method of georeferencing. 
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